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Abstract 

In the last decades, significant progress has been made in understanding the multiphase displacement through porous media with 
homogeneous wettability and its relation to the pore geometry. However, the role of wettability at the scale of the pore remains still 
little understood. In the present study the displacement of immiscible fluids through a two-dimensional porous medium is simulated by 
means of a mesoscopic particle approach. The substrate is described as an assembly of non-overlapping circular disks whose preferential 
wettability is distributed according to prescribed spatial correlations, from pore scale up to domains at system size. We analyze how this 
well-defined heterogeneous wettability affects the flow and try to establish a relationship among wettability-correlations and large-scale 
properties of the multiphase flow. 
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1 Introduction 

The displacement of immiscible fluids through porous media is 
subject of scientific interest, as it is involved in many industrial 
and technological applications, such as oil recovery, water flows, 
and soil treatment. Understanding the dynamics of this process at 
the pore scale is crucial but, on the other hand, very challenging 
due to the many involved parameters, from the competition be- 
tween viscous and capillary forces to the wetting properties and 
porous structure of the substrate. 

In the last decades visualization experiments using simplified 
micromodels of porous media have provided a detailed descrip- 
tion of immiscible flow |U-(5|]. A number of dynamical regimes 
are distinguished in terms of capillary number (ratio of viscous 
to capillary forces), viscosity ratio (ratio of the viscosity of the 
injected fluid to the viscosity of displaced fluid), and wetting con- 
ditions. 

In the limit of high injection rates, the front advance critically 
depends on the viscosity ratio of the fluids. When the injected 
fluid has a lower viscosity than the displaced fluid, the situation 
is highly unstable and ramified viscous fingers are observed. In 
the opposite case, when the injected liquid has a higher viscos- 
ity, the front exhibits a stable displacement with a steady width. 
However, in the limit of very slow injection rates the displacement 
depends entirely on the capillary forces. In this limit, the pressure 
drop due to the viscous flow field can be neglected and the front 
roughens due to capillary pressure fluctuations at the pore level. 
This capillary fingering regime has been successfully described 
by means of statistical models as invasion-percolation an d 
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some more sophisticated pore-network models 1191-11 1H. 

On the other hand, the development of analytical and compu- 
tational techniques for fluid dynamic simulations have provided 
a great advance in the numerical approach to this phenomenon 
1 12W15I1 . However, the continuous description is not sufficient to 
understand the multiphase flow at the pore scale due to the wide 
range of scales involved in this process. To do so, one has to 
comprehend the interfacial dynamics - especially the fluctuations 
generated at the contact line by the substrate heterogeneities- and, 
from this microscopical information, build up a coherent macro- 
scopic description of the flow. 

Mesoscopic simulation techniques provide a robust numeri- 
cal approach with high enough computational efficiency at the 
scale of interest. The fundamental idea is to design simplified ki- 
netic models that coarse-grain the essential microscopical physics 
into a mesoscopic description, so that the averaged properties 
obey the macroscopic transport equations, in particular Navier- 
Stokes and mass diffusion equations. These mesoscale algorithms 
properly reproduce the complex behavior of systems like col- 
loidal suspensions or microemulsions, and have the major advan- 
tage of dealing with changing interfacial topologies in a natural 
way. The most extended ones are direct Monte-Carlo simulations 
(MC) iH], Lattice Boltzmann Method (LBM) [13], Dissipative- 
Particle Dynamics (DPD) 1 1 811 and Multi-Particle Collision mod- 
els (MPC) mm. 

In the present paper we employ the MPC method, originally 
proposed by Malevantes and Kapral [19], to simulate the displace- 
ment of immiscible fluids in a two-dimensional porous medium 
at low capillary numbers. This particle mesoscopic approach is 
an alternative simulation technique that, compared with the most 
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extended LBM, directly includes thermal fluctuations, allows to 
treat systems with an arbitrary number of phases, and provides 
detailed information about the dynamics at the fluid boundaries 
on any scale, as the definition of the physical system is indepen- 
dent on the coarsening process of the mesoscopic approach. 

To deal with multiple phases we adopt the MPC multicolor gen- 
eralization proposed by Inoue et al. 112 111 . We incorporate to the 
model the definition of relative adhesion of the fluid phases to the 
solid phase, in order to control the wetting conditions. Our aim 
is to provide a numerical tool to study the influence of wettability 
on immiscible multiphase flow at the pore scale. Although sig- 
nificant progress has been made toward understanding the effect 
of homogeneous substrate wettability, the role of heterogeneous 
wettability remains still little understood due to the difficulty of 
controlling its spatial definition at the scale of the pores 022l - l25ll . 
Here, we prepare different spatial correlations of local wettabil- 
ity in our numerical simulations and analyze how they influence 
the dynamics, trying to establish a relationship with large-scale 
properties of the multiphase flow. 

In the next Section [2] we provide a description of the MPC 
method including the wettability implementation. In Section[3]we 
analyze the relevant parameters interfacial tension and viscosity, 
and their consistency with thermal fluctuations in terms of cap- 
illary waves. In Section [4] we present our results on multiphase 
flow through porous media analyzing the role of wettability for 
both homogeneous and heterogeneous cases. 

2 MPC model 

The fluid is modeled by a large number N of point-like parti- 
cles, each with the same mass m which move with a continuous 
distribution of velocities. The dynamics of the particles consists 
in a sequence of streaming and collision steps. In the streaming 
step, the particles move deterministically during a time interval 
At according to their individual velocity Vi(t). 

In order to introduce an interaction among particles, i.e. an ex- 
change of linear momentum, they are sorted into collision cells 
£. As we are working in two spatial dimensions, we use a regu- 
lar square grid of typical size a where the number of particles per 
cell, iVf , fluctuates around an average value N c . In every collision 
step, the velocities of the particles are decomposed into the cen- 
ter of mass velocity of the particles in a cell and a remaining, 
fluctuational part. An effective exchange of linear momentum be- 
tween the particles in cell is achieved through a rotation of the 
fluctuational velocity components in each cell. The velocity after 
the effective collision step at time t + At is: 



Vi(t + At) = u^t) + Cl[vi(t) - u s (t)}. 



(1) 



As we assume molecular chaos, the rotation operator £1 in eqn. (Q~|i 
must be statistically independent in each cell and in each time- 
step. Galilean invariance is guaranteed by shifting the coarse- 
graining grid randomly before each collision step [26 1. As linear 
momentum and mass is conserved, the spatially averaged particle 
motion displays hydrodynamic behavior on large length scales. 



Furthermore it has been shown that there is an //-theorem for this 
algorithm 111. 

To deal with multiple immiscible fluid phases we adopt a vari- 
ant of the MPC algorithm proposed by Inoue et al. Kill , that in- 
duces phase segregation. This method has been already success- 
fully applied to study the distribution of droplets in bifurcating 
micro-channels 12711 . In this algorithm, each phase is assigned a 
label (color) c = 1, N v u- At each time step and cell £ the color 
flux 

= - u?) ' Sc,a , (2) 

i=l 

and the color-gradient field 
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are calculated, where V rij is the local gradient of the number 
density n of particles with color j which is estimated from the 
number of particles of the same color in the next-nearest neighbor 
cells. The matrix Kij is the interaction matrix that defines the 
cohesion/adhesion between fluid phases of color i and j. Then 
the multicolor potential energy is defined as 
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and the operator £1 is chosen such to minimize U (£) (in the case of 
2-dimensions this simply means f2(#f) with dU jdO^ = 0, where 
0£ is the angle of rotation respect the center of the cell). Let us 
mention here that this segregation algorithm generates a deple- 
tion layer at the fluid-fluid interface that effectively reduces the 
transport of shear stress through the fluid interface and, in par- 
ticular, a finite slip between the two fluid phases. However, as 
we are mainly working in the regime of small capillary numbers, 
where normal stresses dominate the dynamics at the fluid inter- 
faces, this effect is not relevant. For a more detailed description 
of the model, the reader is referred to [5T|| . 

Thermostatting is required in any non-equilibrium MPC sim- 
ulation due to viscous heating. Here we apply a simple local 
profile-unbiased thermostat that ensures control of the thermal 
fluctuations in the bulk as well as on the solid walls while keep- 
ing unaffected the velocity field ll28ll . The relative velocities to the 
center-of-mass velocity of each cell £ are rescaled after each col- 
lision step as vl = A(£) • uj, where A(£) = ksTN^d/ J^fJi m ' 
(ui — u^) 2 with d dimensions. 

2.1 Wettability implementation 

The multicolor model was originally proposed without speci- 
fying any wetting condition 12111 . However, wettability plays a 
pivotal role in immiscible flow through porous media and, hence, 
needs to be thoroughly implemented in the algorithm. To this end 
we developed a method that accounts for relative adhesion of two 
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fluids to the solid surface while, at the same time fulfilling the 
proper fluid dynamic boundary conditions. 

First, let us review the treatment of the fluid dynamic bound- 
ary conditions for MPC algorithms. In order to guaranty non- 
slip boundary condition for the average fluid velocities, we em- 
ploy t he g eneralization of the bounce-back rule for partially filled 
cells 12911 . In the naive formulation of the bounce back rule, par- 
ticles travel back into the direction of their incidence after having 
collided with the solid boundaries. However, because the posi- 
tion of the solid boundaries relative to the coarse-graining grid 
changes between every interparticle collision step, it is necessary 
to add a virtual phase resting in the walls. This virtual phase takes 
part in the collision procedure to match the bulk particle density 
in the underfilled cells. These wall-particles are generated before 
and removed after every interparticle collision step, and guarantee 
a no-slip boundary condition at the walls II29I1 . 

To define the wetting conditions we include the virtual wall- 
phase in the definition of the interaction matrix, adding an ex- 
tra term that determines the wall interaction for each of the fluid 
phases present in the system. By tuning the relative strength of 
each of them, we can achieve different contact angles of the fluid 
interface at the solid walls. In the following we will refer to a sys- 
tem of two phases, namely fluid (1) and fluid (2), in contact with 
the wall (3), and our matrix is Kij = 1,2,3. 

In order to reduce the slip generated by the phase segregation 
procedure, we set Kis = 1 that fixes the surface tension 013 = 0. 
As there is no distinction between phase (1) and (3) there will be 
no slip for fluid phase (1) at the wall (3). Then, our interaction 
matrix will have the form 



(5) 



According to the definition of the Young's contact angle, re- 
ferred to a droplet of fluid (2) resting on the wall and immersed in 
the ambient fluid (1), we have oyi cos 9 + 023 = C13 that is sim- 
plified to CT12 cos 9 + 023 = with the condition K13 = 1. In the 
remainder of this article we will focus on three particular cases: 
i) For K23 =-lwe find 023 = 012 and therefore a contact angle 
of 9 = 180° while ii) K23 = 1 leads to 023 = and a contact 
angle of 9 = 90°. iii) In the limit K23 > lwe have 023 < and 
therefore 9 ->■ 0°. 

In the following, as we denote the invading phase by (2), we 
will refer to these cases as non- wetting (K23 = — 1). neutral wetta- 
bility («23 = 1) and total wetting, for which we consider K23 = 2. 
The introduction of a negative surface tension for the wetting case 
allows us to avoid the formation of a depletion layer originated by 
the phase segregation mechanism and avoid any slippage of the 
invading phase (2) at the wall. In Fig. Q] we show the equilibrium 
configurations of a droplet of phase (2) between two solid parallel 
walls for the three wetting conditions we consider. 
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Figure 1: Equilibrium contact angle for the three considered wetting 
conditions in a planar geometry (from left to right): non-wetting, neutral 
wettability and total wetting. 



3 System parameters and thermal fluctuations 

In this section we analyze how macroscopic quantities, like the 
dynamic shear viscosity and interfacial tension, emerge from the 
definition of our microscopic system parameters, that are the cell 
size a, time step At, particle mass m, mean number density N c 
defined as the average number of particles per cell, temperature 
temp (in UbT units), and the entries Kij of the interaction matrix. 
We normalize At = 1, a = 1, m = 1, and set N c = 10. For 
all the simulations we fix temp = 1CP 3 , which implies a mean 
free path for a particle in between two consecutive collisions A = 
At^temp/m « 0.03. 

The behavior of a single phase does not depend on the strength 
of the self-interactions among particles of the same species, so 
we can normalize the diagonal terms of the interaction matrix to 
Ku = 1 without loss of generality. Fitting the velocity profile 
of a simple fluid in a planar channel to the Poiseuille flow we 
estimate the value of the dynamic viscosity of the fluid phase to 
be 11 = 0.750 m/aAt for the set of parameters considered. 

Since the interfacial tension the fluid phases (1) and (2) is not a 
direct input for our algorithm but is controlled by the interaction 
term K12, we estimate it by terms of the Laplace law. Initially 
a square droplet of fluid (1) is placed in the center of a simula- 
tion box filled with fluid (2). After a transient time, it relaxes to 
a configuration fluctuating around a circular shape. We measure 
the radius R of the averaged droplet contour and the correspond- 
ing difference AP = P\ — P2 between the pressure in the droplet 
(1) and the in the ambient fluid (2) for different values of the in- 
teraction term K12. The pressure is computed by means of the 
microscopic definition of the local stress tensor 



rrii ■ Vi a ■ Vip 



VAt ^ 

i—l 



Ap ia ■ r if s , (6) 



where Vi is the velocity of the particle before the collision step, 
Ap/mi = Vi(t+At)~Vi(t) and r*i is the particle position referred 
to the center of the cell. The expression eqn. (|6) for the stress 
tensor is given, e.g., in Ref. i30ll . 

Plotting AP against the radius of the averaged droplet contour 
in Fig. |2] shows the expected linear relation AP oc 1/P accord- 
ing to Laplace's law. The proportionality constant is the value 
of the interfacial tension a. We observe that the interfacial ten- 
sion becomes non-zero as soon as H12 ^ K\i = ^22- After a 
short rise while decreasing K12 , it describes a plateau for values 



3 



K12 < — 1. According to this, in forthcoming simulations we fix 
our fluid-fluid interaction K12 = — 1, for which we have an inter- 
facial tension a = (0.0011 ± 0.0004) m/At 

Finally, we check that the thermal fluctuations are consistent 
with the estimated value of the interfacial tension analyzing the 
capillary waves observed on a planar fluid-fluid interface. An ini- 
tially straight interface will roughen by the motion of thermally 
activated capillary waves. In the capillary wave spectrum, each 
Fourier component hk = Y^rjZo hj{t)e 2m i k / L of the interface 
displacement contributes according to the equipartition theorem 
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(7) 



leading to an interface roughness proportional to ^ksT/a yl] 



3211 . In Fig.|2]we plot the average amplitudes (\hk\ 2 ) of the planar 



capillary waves as function of the wave number k. As expected 
from the assumption of a wave number independent of the inter- 
facial tension a we observe that the spectra for different system 
sizes collapse for the value a = (0.0023 ± 0.0008)m/At, which 
is in good agreement with the estimate obtained from the Laplace 
law for circular droplets of different radii. Therefore, the model 
displays the correct thermodynamic behavior and interfacial fluc- 
tuations. 



4 Flow in porous media 

Immiscible displacements in porous media with both capillary 
and viscous effects can be characterized by two dimensionless 
numbers, the ratio M = [i& of the dynamic viscosities \i of 
the invading fluid (i) and the displaced fluid (d), and the imposed 
capillary number, i.e., the ratio of viscous to capillary forces given 
by Ca = av/fii, where a is the interfacial tension of the fluid 
interface and v is the average front velocity. 

When viscous forces dominate (Ca > 1), the invading front ex- 
hibits either stable advance or viscous fingering (M > 1 and M < 
1 respectively). However, at low Ca the displacement is solely 
governed by capillary forces This propagation mechanism 

can be captured using concepts from invasion-percolation JfjrU I- 
In the following, we consider a porous matrix which is completely 
filled with the initial fluid phase (d). By means of an applied pres- 
sure gradient a secondary fluid phase (i) is pushed into the matrix 
displacing phase (d). 



4. 1 System setup 

The porous matrix is represented by an assembly of non- 
overlapping circular disks with fixed position whose preferential 
wettability by the invading fluid (i) can be selected from the three 
different cases: non-wetting, neutral and wetting. Our porous re- 
gion is a two-dimensional rectangular channel with length L and 
width H that extents < x < L in the direction of the applied 
pressure gradient. Periodic boundary conditions are considered 
into orthogonal direction. 

To generate the pressure gradient that drives the flow we de- 
fine an inflow region x G [— /, 0] being free of obstacles where 
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Figure 2: Interfacial tension analysis. Top: Laplace AP = P\ — P2 
plotted vs. the interfacial curvature 1/R for droplets of different size. 
In the inset we plot the interfacial tension o estimated from the slope 
of the pressure/curvature plot vs. K12. Bottom: Power spectrum of 
thermally excited capillary wave amplitude for a planar interface. The 
spectra for different system sizes collapse according to eqn. for 
o = (0.0023 ± 0.0008) m/At, in agreement with the estimate from 
the Laplace pressure. 



a constant body force g acts into the x direction on the particles. 
The extension of this inflow region is / w 0.3L into the direction 
of the flow. As a result of the body force in the inflow region a 
pressure gradient builds up that drives the fluid particles through 
the obstacles. 

The value of g is constrained to be small enough in order to 
avoid compression of the fluids, i.e. a Mach number Ma = v/c < 
1, where c is the speed of sound. In the following we will apply 
g = O(10~ 5 ), which gives us Ma = 0(1O~ 2 ), Reynolds num- 
bers Re = pvC/fi = O(10 _1 ) with C the typical distance be- 
tween the obstacles and mass density p, and Ca = O(10 _1 ) for 
the system sizes and the set of parameters we consider, as intro- 
duced in section [3] 
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4.2 Homogeneous wettability 

First we analyze the dynamics of the forced displacement in 
a substrate composed of monodisperse disks. Previous experi- 
ments in quasi-2d micromodels ||51lZ|], as well as numerical ap- 
proaches 



3311 clearly indicate that wetting conditions strongly 



influence the displacement of immiscible fluids through porous 
media in the capillary fingering regime. 

These analyses also evidence that, when viscous forces are not 
strictly zero, the flow may be influenced by the pore geometry. 
For instance, in the case of a broad throat-size distribution the 
viscous pressure drop may become on the order of magnitude of 
the difference of Laplace pressure between a small and a large 
throat, what makes the invasion-percolation description not longer 
valid 
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To discern the effect of wettability from effects related to pore 
geometry and from dynamic front instabilities (i.e., the Saffman- 
Taylor instability) we strictly consider equal dynamic viscosities 
M = 1 and porous matrices with low packing density (j> = A a /A, 
where A Q is the area occupied by the obstacles an A is the total 
area of the simulation box. In the following we employ a ran- 
dom configuration of circular obstacles with radius r, generated 
by means of random sequential addition 83411 . for which we im- 
pose = 0.15 and a minimal separation distance d = 2r between 
the obstacles. 

Figure [3] displays snapshots of the invading front for the three 
different cases of wettability, as well as details of the flow around 
one single obstacle for each case. Initially the matrix is com- 
pletely filled with phase (d) and then the invading fluid (i) is driven 
from the left by means of the pressure gradient generated in the 
inflow region. The applied body force on particles within this 
region is g = 1.0 x 10~ 5 , what generates a pressure gradient 
AP « 0.008 m/At 2 . In these figures we show snapshots of the 



invading front close to percolation. They correspond to different 
simulation times as the velocity of invasion clearly depends on 
the wettability. 

In the non-wetting case we observe that, almost from the very 
beginning, the front roughens and fingers are formed. The fluid 
in these fingers pinches off leaving behind a number of discon- 
nected domains of the displaced fluid which are trapped between 
the disks. This gives rise to a considerably high residual satu- 
ration as we observe in Fig. |4] where the evolution of saturation 
S = A/Aq is depicted versus time, with the area A of displaced 
phase (d) and total area Aq of the void space. On the other hand, 
the formation of new fluid-fluid interface slows down the dynam- 
ics, as we observe when plotting the flux measured at the end of 
the porous matrix for both fluid phases. The flux of the invading 
fluid, instead of remaining constant in time, as would correspond 
to a stable front, clearly decreases in time as the invading fluid de- 
scribes intricate paths through the porous structure. These plots 
also shows that even after breakthrough, when the flux of the dis- 
placed fluid displays a first peak, some other fingers still develop, 
and we observe a secondary percolation of the invading fluid be- 
fore reaching the steady state at which the flux becomes constant. 

In the case of neutral wettability the dynamics of the invading 
front is considerably different. We observe much smaller trapped 
clusters of the displaced fluid. A similar behavior has been re- 
ported in Ref ^ for quasi-2d experiments. This is clearly mani- 
fested by the lower value of residual saturation in Fig. |4] However, 
in this graph we observe that the velocity of invasion, that we de- 
fine as the velocity of an equivalent stable front with v ~ LdS/ dt 
for the same imposed Ca, is very similar for both the non-wetting 
and the neutral case until the curves saturate. As in both situa- 
tions the work injected into the system can be invested either in 
viscous dissipation or in formation of new fluid-fluid interfaces, 
this means that the length of the fluid-fluid interface should be 
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also comparable in both cases during this time span. To figure 
this out and complement this analysis, we depict in Fig. [5] how 
the length of the fluid-fluid interface evolves in time. Certainly, 
we observe that the length growth of the fluid-fluid interface is 
very similarly for both cases during this time span, after which 
it practically remains constant for the non-wetting case while is 
drastically reduced, practically constant, for the neutral case. 

Finally, we consider the case of perfectly wetting obstacles for 
the same imposed pressure gradient as in the previous cases. In 
this case the front advance is clearly stable, and the flux of the 
invading phase fluctuates around a constant value as well as the 
length of the fluid-fluid interface remains constant until break- 
through. In this case the time until percolation is much shorter, 
which is readily explained by the gain in wetting energy in addi- 
tion to the pressure gradient. After percolation, virtually all the 
initial saturation is taken out as can be seen in Fig. [4] 
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Figure 4: Saturation of the displaced phase vs. time for homogeneous 
wettability: non-wetting (red), neutral (green) and wetting (blue). Re- 
sults are shown for low driving g = 1.0 x 10~ 5 (full symbols) and higher 
driving g = 1.5 x 10 (open symbols). The capillary number for each 
case is specified in the inset. Dashed line is depicted for comparison with 
Fig.\5\ Top: Evolution of the flux of the invading (black) and displaced 
fluid phase ( red) during the invasion for the case of non-wetting (left) and 
wetting (right) pore walls for the low driving case. The dashed and dotted 
lines indicate breakthrough and secondary percolation, respectively. 



Next we analyze how the front advance depends on the applied 
driving g. In accordance with experimental results Jfi we observe 
shorter percolation times and lower values of the residual satura- 
tion for higher driving g = 1.5 x 1CP 5 . For comparison with the 
previous results at lower Ca, we include the evolution of S vs. 
time in Fig. |4] It is clear that, in the case of higher Ca, the capil- 
lary fingers can overcome more easily the local Laplace pressure 
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Figure 5: Evolution of the fluid-fluid interfacial length in time for 
the three different homogeneous wetting conditions considered: non- 
wetting, neutral and wetting from top to bottom. In the insets we show 
a detail of the fluid-fluid interface for the non-wetting (top) and wet- 
ting (bottom) cases, corresponding to the snapshots presented in Fig\3\ 
Dashed line is depicted for comparison with Fig. [4] 



required to invade the pores, that in this regime is basically given 
by the pore geometry (at least in the limit M = 1). This clearly 
favors the extraction of initial fluid, as we can appreciate in Fig.|H 
where the residual values are smaller than the ones obtained for 
any case of wettability at smaller Ca. 

4.3 Heterogeneous wettability 

In this section we analyze how wettability heterogeneities at 
pore scale affect the dynamics of fluid invasion. Previous works 
dealing with heterogeneous wetting conditions are present in the 
Literature, for instance experiments dealing with a mixture of 
wetting and non-wetting beads 13511 or composite porous medium 
made of blocks of the same permeability but opposite wettabil- 
ity II25II . We find also detailed numerical approaches, like simula- 
tions of immiscible displacement through a realistic pore structure 
whose wettability distribution is taken from microtomographic 
images of reservoir rocks |]24|. 

However, our aim here is to analyze the effect of such hetero- 
geneities using well-characterized spatial distribution of wettabil- 
ity down to the pore scale. This implies being able to indepen- 
dently assign to each pore a well defined wettability according to 
a predefined spatial distribution. 

In order to identify the pores in our porous matrices we con- 
sider a Delaunay triangulation based on the centers of the ob- 
stacles. According to the construction no obstacle-center is in- 
side the circumcircle of any other triangular pore. Once defined 
the pores, we independently assign wetting conditions to each of 
them. This means that our obstacles are divided in sectors accord- 
ing to the pore definition, and different wettabilities are assigned 
to each of them, to build up an heterogeneous pattern at the pore 
scale. 
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Figure [6] shows the distribution of two immiscible fluids of an 
initially random mixture in presence of a pore-scale wettability 
pattern after spontaneous phase separation. It can be observed 
that each of the two immiscible fluids is in contact solely to the 
highly wettable parts of the pore wall. Deviations from an ideal 
distribution can be easily explained by the presence of thermal 
fluctuations. 




Figure 6: Pores defined in terms of Delaunay triangularization and wet- 
tability definition at pore scale. On the right we show the final state 
of a 1:1 binary mixture of two immiscible fluids after relaxation, which 
matches with the imposed wettability pattern defined in terms of wetting 
and non-wetting pores. 

Next we consider the same porous geometry than that of sec- 
tion l4.2l and impose wetting patterns with well defined spatial cor- 
relations, in terms of wetting and non-wetting domains with typ- 
ical size £. To generate these patterns we start from a randomly 
distributed wettability for each pore and introduce correl ations by 
means of a simple stochastic procedure in which 

- we randomly choose a pair of pores 

- we assign to the second pore the wettability of the first one 
with probability p(£) that depends on the distance £ between 
both of them, i.e., the distance between the incenters of the 
triangles. 

We choose this function to be the error function p(£) = 
3 fl + erf (^7§j)> that exhibits a sharp transition between 1 
and at a typical distance /i. We fix a = 3yU. Due to compu- 
tational system-size limitations, the number of pores in our sys- 
tem is not large enough and there is some deviation between the 
typical distance ii imposed in the pattern generation and the spa- 
tial correlation length £ that the patterns finally exhibit. However, 
when fitting the spatial correlation of the samples to the functional 
form C(£) ^ exp(—£/£ l ) we obtain a correlation length £ that de- 
pends linearly on the \i parameter as £ = 1.2 + /j, (goodness of fit 
C = 0.98). On the other hand, the ratio the total area occupied 
by the wetting and non-wetting pores, A w and A nw respectively, 
is approximately constant for all the patterns A w /A n - W ss 0.5. 

In the following we will analyze how the variation of the spatial 
correlation of our wettability patterns affect the invasion when 
considering the same arrangement of obstacles than that of the 
homogeneous case. 



In the case of spatially uncorrelated wettability, i.e. small val- 
ues of Li, the front needs to invade many small non-wetting do- 
mains in order to percolate. We show these pore -breakthrough 
events in Fig. [8] for a given realization at small domain size with 
ll = 4. For the sake of better visibility the invaded non-wetting 
pores are marked in red. 

If we increase the correlations, percolation of the invading fluid 
requires filling of larger non-wetting domains. This slows down 
the invading process as can be seen in Fig. [7] where we show 
the evolution of saturation S for different values of the ii param- 
eter. On the other hand, invading these larger non-wetting do- 
mains is more expensive in terms of surface energy, and therefore 
many of the non-wetting domains remain unexplored. As we see 
from the plots shown in Fig. [7] the residual saturation increases 
monotonously with the correlation length for this range of inter- 
mediate correlations. 

Finally, for even larger correlation lengths we observe that only 
few bottle-neck non-wetting pores need to be invaded to achieve 
percolation and the invading fluid rapidly floods the large wetting 
domains. Now, the residual saturation is virtually given by the 
ratio A w /A n _ w and does not depend any more on the correlation 
length. 

We have to mention here that these results are statistically poor 
due to the computational limitations to go for larger systems and 
have better defined spatial correlations. However, they are still 
good enough to show qualitatively that the fluid invasion is sensi- 
tive to the spatial distribution of wettability at pore scale. In Fig. [7] 
we show the dependence of the residual saturation after the inva- 
sion on the correlation parameter fj,, where each value has been 
averaged over several equivalent samples and independent real- 
izations on each of them. The increase of the residual saturation 
with the correlations of the wettability patterns is clearly appre- 
ciable, as well as the subsequent plateau due to the presence of 
bottle-neck non-wetting pores. 

Finally, we conclude this analysis increasing the imposed Ca, 
as we did in the case of homogeneous wettability. Now the higher 
pressure gradient washes out the effect of the wettability patterns 
on the flow advance. We appreciate in Fig. [7J that the saturation 
curves for different /i are similar and the difference in residual 
saturation diminishes, in comparison with the lower driving case. 

5 Conclusions 

We have adapted the MPC multicolor algorithm to control the 
wetting conditions at the solid boundaries in order to simulate 
immiscible displacement trough porous media. This is an alterna- 
tive particle-based simulation approach that provides microscop- 
ical details of the dynamics at the fluid boundaries and directly 
incorporates thermal fluctuations. 

As a first application of this method, we analyze the depen- 
dence of the multiphase dynamics on the spatial correlations of 
wettability, analyzing the evolution of the initial saturation and 
particle flux for different wetting conditions. 
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Figure 7: Saturation as function of time for heterogeneous wettability 
patterns with different /i (single realization for each case). Top: At inter- 
mediate correlation lengths, we observe a longer time until percolation 
accompanied by a higher residual saturation for increasing £. The curves 
correspond to parameter /j, = 4, 5, 6, 7 (from bottom to top). In the in- 
set the saturation curves for larger correlation of the pattern are depicted 
(p, — 8, 9, 12). In this case the invading front rapidly percolates with con- 
stant value of residual saturation. All the simulations are executed with 
driving g — 1.0 x 10 . Bottom: Saturation curves for higher imposed 
Ca (g = 1.5 x 10~ 5 ), for the same parameters values fj, = 4, 5, 6, 7. The 
higher applied pressure gradient washes out the effect of the wettability 
patterns, and both the invading velocity and the residual saturation be- 
come similar for all the cases. Inset: Residual saturation vs. /i averaged 
over several equivalent patterns for each value of /i and independent re- 
alizations on each of them. Circles and squares correspond to the lower 
and higher imposed Ca, respectively. The higher driving clearly dimin- 
ishes the effect of the wettability pattern. Lines are included to guide the 
eye: y ~ 0.03/1 (solid) andy ~ const, (dashed) 



different prescribed spatial correlations at the pore scale. We ob- 
serve that the dynamics is influenced by such spatial correlations. 
While in the range of intermediate correlations the residual sat- 
uration increases with £ whereas the invading process becomes 
slower, for large enough correlations (of the order of the system 
size) the presence of some bottle-neck pores controls the dynam- 
ics, and the invading fluid rapidly percolates, while the residual 
saturation remains constant given by the ratio of wetting and non- 
wetting pore area. 

These qualitative results manifest the interest of a further explo- 
ration of the influence of the wettability distribution at the scale 
of the pore. In this work, the statistics of our heterogeneous pat- 
terns is limited by the system size, that is in turn limited by the 
computational costs of the MPC algorithm implemented in se- 
quential mode. However,this algorithm provides a good level of 
parallelism. In the free streaming step, the particles are prop- 
agated according to their velocity without interacting with each 
other. On the other hand, the multi-particle collision step is per- 
formed cellwise. This means that the necessary calculations are 
completely independent from cell to cell, so that they can be also 
executed in parallel. This will be certainly exploited in further 
works in order to carry out more accurate simulations to analyze 
microscopic flow mechanisms that are observed at pore scale. On 
the other hand, this MPC multicolor model has the advantage of 
being able to deal with an arbitrary number of phases, what is a 
interesting feature to simulate more realistic systems. 

We thank Yasuhiro Inoue and Anikka Schiller for helpful dis- 
cussions as well as Stephan Herminghaus for his comments and 
support. Funding from BP Exploration Operating Company Ltd. 
within the ExploRe research program is gratefully acknowledged. 
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